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CHAPTER 1 


INTRODUCTION 

There is a current emphasis in computational fluid dynamics (CFD) to develop tech- 
niques for solving flow fields about geometrically complex bodies in two and three dimen- 
sions. A major difficulty in the computation of such flows is the specification of a proper 
discretization of the flow field about the geometry of interest. A numerical solution of the 
governing partial differential equations for the flow about a given body requires the gen- 
eration of an appropriate grid (a set of points or volumes) spanning the region of interest. 
The differential equations are then approximated by a set of coupled, algebraic, difference 
equations, and the system is solved for functional values at each point or volume. Tradi- 
tionally, a single curvilinear grid is generated for the entire flow-field region, and is often 
“body-oriented” (body surface corresponds to a constant grid-coordinate surface) to aid in 
the application of boundary conditions and/or to permit a simplification of the governing 
equations. Problems with complex geometries, however, pose difficult challenges for grid 
generators; the creation of a single grid of acceptable quality about the entire configura- 
tion may even be impossible. In fact, Chapman (1980) and Kutler (1985) have identified 
this problem of three-dimensional grid generation as a primary pacing item for the field of 
CFD. 

The discretization of a flow field directly affects both the smoothness and the accuracy 
of a numerical solution; hence, attention must be given to the quality of the chosen grid. 
Unfortunately, there are few quantitative measures of the “goodness” of a mesh and most 
assessment, in practice, is done by eye. Nevertheless, there are some important, qualitative 
considerations. 

To best represent the flow-field solution, the grid should be clustered in expected 
high solution-gradient regions. However, the efficient use of computer resources precludes 
the employment of a fine mesh throughout the field but suggests the use of nonuniform 
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grids with gradual stretching from coarse areas to the desired fine areas. Thompson et 
al. (1982) show that the truncation error in difference expressions is a function of local 
grid spacing and rate of change of the spacing, as well as the departure of the grid from 
orthogonality. To obtain a numerical solution, derivatives along the mesh coordinate lines 
must be evaluated. Therefore, abrupt changes in either grid-line orientation or spacing will 
lead to inaccuracies in the estimation of these derivatives. Often these inaccuracies are 
manifest as local “kinks” in the solution contours. Smoothly varying grid coordinate lines 
and spacings are highly desirable, as are nearly orthogonal coordinate-line intersections. 
Furthermore, Steger (1982) warns that grid properties such as smoothness, skewness, and 
mesh cell aspect ratio may not only affect the accuracy of the solution but also the con- 
vergence rate of the numerical algorithm. Some algorithms appear to be very sensitive to 
undesirable mesh properties (e.g., successive line overrelaxation iterative methods and cell 
aspect ratio variations; the Beam- Warming algorithm and highly skewed grids). 

The foregoing considerations impose stringent constraints on the task of grid genera- 
tion; consequently, the development of sophisticated algebraic and differential techniques 
for grid generation is an active area of research. But multiple-body and geometrically com- 
plex problems, especially in three dimensions, are difficult with even the most advanced 
methods. Often the resulting grid about the configuration (after significant investment in 
research time) is of compromising quality. 

A perhaps obvious alternative to the specification of a single grid in a complex region 
is the use of some sort of blocked or zonal approach to grid generation. The zonal concept 
is intended to simplify mesh generation by a subdividing the flow field and generating 
grids within each subregion. Researchers are exploring block approaches for three primary 
reasons: (l) to permit the solution of problems which are intractable with a single-grid 
approach because of their geometric complexity; (2) to increase solution accuracy by se- 
lectively refining subregions of a previously specified, global, single grid; (3) to facilitate 
a block-processing computational approach wherein information for only one zone needs 
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to reside in the computer core at any given time, thereby relaxing memory limitations, 
or where computers operate in parallel to compute the flow fields in various zones. The 
use of such a zonal approach requires the development of a technique for treating the grid 
interface boundaries. Ideally, the procedure used to transmit information from subregion 
to subregion will maintain a stable, accurate, and conservative computation. 

There are two general classes of blocked approaches to grid generation and solution 
methodology: overlapping grid systems and patched grid systems (point-continuous and 
point-discontinuous). Figure 1 illustrates these concepts for a simple, two-dimensional, 
single-body problem. The salient features of each class are discussed in the following 
paragraphs. 

The overlapping approach (fig. la) provides the most flexibility in grid generation of 
the blocked schemes. A mesh is generated about each body or component of a complex 
configuration with the sole constraint that there be some degree of overlap between neigh- 
boring grids. This method is especially attractive for time-dependent applications in which 
one body is in unspecified motion with respect to another (such as in a store-separation 
computation). Atta (1981) and Atta and Vadyak (1982) test this approach with the full 
potential equations in two and three dimensions for steady-state problems. They employ 
either Dirichlet or Neumann boundary conditions at mesh interfaces using interpolated 
values from neighboring grids. Similarly, Thompson (1981) uses an overlapping grid sys- 
tem in computing transonic wing/body/store flow fields by solving the small disturbance 
potential equation. More recently, Steger et al. (1983), Benek et al. (1983), Benek et al. 
(1985), and Dougherty (1985) carefully outline the data structure details for an overset 
grid scheme they refer to as ‘‘chimera.” It is clear from their work that the overlapping 
approach greatly simplifies the task of grid generation but does complicate data manage- 
ment. Also, Benek et al. (1983) document problems with the scheme as a shock wave 
passes through the grid boundaries and speculate that the cause is the lack of conservation 
in these regions because of the nonconservative interpolations at the interfaces. There is 
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further evidence that this may indeed be the case based on the work of Hessenius and Pul- 
liam (1982) and Rai et al. (1984). Both references cite examples of computations in which 
the introduction of a degree of nonconservation in the vicinity of a shock wave produces a 
converged but incorrect solution (the shock does not assume the proper strength, shape, 
and/or position). Berger (1984) outlines a conservative treatment of arbitrarily oriented 
overlap boundaries for two-dimensional computations by performing a flux balance at the 
interface. Unfortunately, the generality and flexibility of the overlapping scheme geometri- 
cally complicate the procedure of balancing fluxes. To date, Berger’s method has not been 
implemented. 

Overlapping grid systems, while providing a simplified treatment of complex geome- 
tries, are used for other reasons as well. Berger (1982) develops a technique for hyperbolic 
problems for automatically inserting or deleting overlapping grids based on local truncation 
error, thereby increasing the accuracy of the computations in an efficient manner. Caruso 
(1985) modifies Berger’s approach for use in solving elliptic flow problems. Similarly, Lom- 
bard and Venkatapathy (1985) examine the use of overset grids for improved resolution of 
flow field features. Chang et al. (1985) overcome limitations in computer core memory by 
zoning the internal flow field in the Space Shuttle main engine using point-continuous but 
overlapping meshes. 

There are two main variations of the patched grid approach to numerical compu- 
tations: point-continuous and point-discontinuous (see figs, lb and lc). Both types of 
patched grid systems permit a segmentation of the flow field into zones, but require that 
neighboring zones meet along a common boundary (curve in two dimensions or surface 
in three dimensions). The point-continuous system, however, also requires continuity of 
grid lines across the patch boundaries. Two-dimensional examples of the point-continuous 
solution approach are found in the work of Lasinski et al. (1982) in computing the flow 
about a tri-element airfoil and in the work of Norton et al. (1983) in calculation of the 
flow through turbine cascades. Rubbert and Lee (1982) use a point-continuous patching 
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approach for the generation of three-dimensional grids about aircraft configurations, as do 
Weatherill and Forsey (1984) for wing-canard computations. 

The point-continuous patched schemes referenced above, while permitting compo- 
nent adaptive grids for each body or component part, still necessitate the development 
of fairly sophisticated generation techniques because of the enforcement of pointwise con- 
tinuity at the interfaces. (The grids within each block cannot be generated completely 
independently.) Furthermore, the pointwise-continuous constraint invariably forces the in- 
troduction of one or more mesh singularities which require special treatment in the flow 
solver. However, the “interpolation-free” zonal interface conditions make the computation 
naturally conservative (an attractive feature), unlike other methods discussed thus far. 

The point-discontinuous approach has been used in a grid refinement mode by Berger 
and Jameson (1984), Baker et al. (1985), Holst et al. (1985), and Eriksson (1985) wherein a 
global coarse mesh is selectively refined in “high action” regions (usually by the subdivision 
of a single coarse grid cell into several fine grid cells). However, it may still be a difficult and 
time-consuming task to generate the initial coarse grid. Nevertheless, localized patching 
provides the desired resolution in important areas of the flow. 

A more general point-discontinuous, patched grid scheme (fig. lc) allows for indepen- 
dent discretization of each subregion or zone; often undesirable mesh singularities may be 
completely avoided. If the zones are judiciously chosen to be rather simple regions, existing 
grid generation techniques (often algebraic) can be used. As with overlapping grid systems, 
special boundary procedures must be developed to handle the mismatch of grid lines at 
the zonal interfaces. Cambier et al. (1984) and Veuillot and Cambier (1984) introduce a 
characteristic boundary method for patching two-dimensional, point-discontinuous grids. 
Because their procedure is nonconservative in nature, they have fit the shock waves in tran- 
sonic flow computations. The resulting solutions are smooth and continuous across the 
patch boundaries. Bush (1985) also uses a characteristic (and therefore nonconservative) 
patching method in two dimensions to compute the flow about an external compression 
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inlet. 

Rai (1986) outlines a fully conservative treatment of zonal boundaries for two dimen- 
sional Euler equation computations on discontinuous, patched grid systems. Success with 
this method has been clearly demonstrated in the solution of both steady-state and time- 
dependent problems (two-dimensional rotor-stator interactions) in subsonic and supersonic 
flows (Hessenius and Rai, 1986; Rai, 1985a; Rai, 1985b; Rai, 1985c). Discontinuities move 
freely through the grid interfaces to assume their proper locations, a consequence of the 
conservative nature of the scheme, and the solution is continuous across these boundaries. 

Rai’s approach to patched grid computations has the desired qualities of a zonal 
scheme. By subdividing the flow field into well-shaped regions or zones, grid generation for 
complex geometries is greatly simplified and may be accomplished with existing techniques. 
The zones permit a block-processing approach to computation. The interfacing procedure 
is shown to be stable and accurate and is fully conservative in the discrete sense. 

The purposes of this work are (1) to establish the usefulness of Rai’s (1986) technique 
in the solution of geometrically complex problems in two dimensions, and (2) to extend the 
zonal scheme for use in three dimensions and demonstrate the procedure with a multiple- 
body application. To this end, Chapter 2 contains a review of the two-dimensional method 
as outlined by Rai (1986), while Chapter 3 presents computational results for flow about 
a double-airfoil configuration using two different patched grid systems. The development 
of the three-dimensional technique comprises Chapter 4. Demonstration computations of 
flow about a wing-canard combination are contained in Chapter 5. Concluding remarks 
and recommendations for further study are given in Chapter 6. 
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Fig. 1 




(c) 

ck approaches to grid generation (a) overlapping 
patched: point-continuous (c) patched: point-discontinuous 
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CHAPTER 2 


REVIEW OF RAI’S TWO-DIMENSIONAL 
PATCHED GRID SCHEME 

Rai’s (1986) technique permits the patching together of subgrids to form a global, two- 
dimensional grid (see fig. lc). Neighboring subgrids share a common “zonal” boundary (a 
coordinate line) , but along this boundary the grids may be point discontinuous. The proper 
treatment of the zonal-boundary points is crucial to the maintenance of a conservative, 
accurate, and stable computation. Rai’s method ensures global conservation of mass, 
momentum, and energy in a discrete sense at the zonal interface. (A scheme is globally 
conservative if a summation of the fluxes throughout the flow field results in a sum of 
only the boundary fluxes (i.e., all interior fluxes cancel) with no residual fluxes at a zonal 
interface.) 

This technique (Rai, 1986) is developed within an explicit, first-order accurate (in 
time and space), finite-difference framework for the Euler equations and will be outlined 
as such in this chapter. However, the method has been extended for use with second-order 
accurate, implicit integration schemes (Rai, 1985). 

Consider the patched grid of figure 2. This two-zone grid, sketched for a supersonic 
blunt-body application, contains a single zonal boundary along a constant ^-coordinate 
line. Note that across this line (line AB) the grids are discontinuous and created indepen- 
dently. 

The method of solution on this patched grid may be considered a three-step process: 

1) The interior grid points in zones 1 and 2 are updated by integration of the governing 
equations in any conventional manner. 

2) The zonal-boundary points in one of the two zones are updated using a procedure 
that guarantees conservation across the boundary. (The choice of the zone to update 
first does not affect the conservation properties of this method.) For example, if one 
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chooses to update the zone-2 points on the zonal boundary, the procedure to update 
point (j,l) (see fig. 2) involves an integration of the governing equations using the 
flux cell RSTU. Fluxes at sides UR, RS, and ST are available from zone 2. The flux 
at side UT is determined by an integration of the zone-1 flux distribution along the 
line CD between points U and T. 

3) The solution at points on the zonal boundary in the neighboring zone (in this case, 
zone 1) is found by interpolating the solution from step 2 along this boundary. 

2.1 Solution at Interior Points 

The unsteady Euler equations in two dimensions may be written as 

Qt + E x + F y = 0 ( 2 . 1 ) 


where 


Q = 


p \ 

pu 

E = 
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and 


P = (7- 1)[«- 2 P( u2 + t;2 )l 


(Note that p is the density; p the pressure; u and v are the velocity components in the x 
and y directions, respectively; and e is the total energy per unit volume.) Applying the 
independent-variable transformation 

r(*> = t 


= « <0 (*,»,«) 

>;<*> = i 



for zone 1 
for zone 2 


to equation (2.1) for each zone of figure 2, one obtains the Euler equations in generalized 
coordinates 


where 
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A conservative finite-difference scheme for integrating equation (2.2) is given by 
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= 0 


(2.3) 


where F, + 1 / 2 ,/t and Fy^ 1/2 are numerical fluxes consistent with the transformed fluxes 
E and F. An evaluation of the numerical fluxes F and F at time-level n in updating 
to Q n+1 would characterize an explicit finite-difference scheme. The following discussion 
assumes the use of such a scheme. 

The interior points of both zones 1 and 2 (fig. 2) are explicitly updated using equation 
(2.3) with centrally differenced metric terms. Points along the zonal boundary are updated 
as outlined in the following section. 


2.2 Solution at the Zonal Boundary 

Consider the point (j , 1 ) of zone 2 on the zonal interface and its associated cell area, 
RSTU (fig. 2). (Note that the fluxes in eq. (2.3) are defined at half points and are, 
therefore, evaluated at cell sides SR, ST, TU, and RU.) This cell area is found by an 
extrapolation of the coordinate lines of region 2 to a previously determined flux-balance 
line (labeled CD) in the neighboring zone 1. The flux balance line is chosen halfway 
between the last two constant ^-coordinate lines in zone 1. The metrics at the point (j,l) 
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are then defined 


= 2^(*j+i,i ~ 

i.i = 2^(»y+i,t “ yy-1,1 

i x v)j, i = ^( x >,3/2 — x i > 1 / 2 ) 

~ ^(yj.3/2 ~ t/j, 1 / 2 ) 

Metrics at cell sides (half points) are obtained by averaging the above expressions between 
neighboring grid points. 

Three of the required four fluxes necessary to update the point (j,l) of zone 2 us- 
ing equation (2.3) are constructed using only zone-2 information. The ^-direction fluxes 
{Ej- x/2,i an d Ej+ 1/2,1 through cell sides RU and ST, respectively) are formed using the 
metrics x ^ and y, at the half points and values of the dependent variables at points (j- 1,1), 
(j,l), and (j+1,1). The 77-direction flux F J)3 / 2 through cell side SR is computed using the 
^-difference metrics (at half points) and values of the dependent variables at points (j,l) 
and (j,2). Thus far, this is the standard procedure for updating interior points as well. 
However, the flux into the region 2 cell at point (j,l) from region 1 (the Fj,i/2 term) must 
be conservatively interpolated from known zone-1 information in the following manner. 

First, one constructs the numerical flux in zone 1 at all points along the flux balance 
line using values of the dependent variables in zone 1 along the last two constant 77- 
coordinate lines and zone-1 metrics. These fluxes are then divided by an appropriate unit 
of length, the metric corresponding to a zone-1 cell side along the balance line, to create 
flux per unit length. The unit of length for zone-1 point l is best thought of and computed 
as the arc length along the balance line between 1 + 1/2 and l — 1/2. Then assuming a 
piecewise-constant variation of the zone-1 numerical flux along the flux balance line (see 
fig. 3; l points (•) are from zone 1 and j points (x) are from zone 2), one obtains the 
flux into the zone-2 zonal-boundary cell by an integration of the zone-1 flux distribution 
along the appropriate portion of the flux-balance line. (Consider s a running parameter 
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along this line: increasing s with increasing £.) Therefore, the flux from zone 1 into the 
cell RSTU (fig. 2) is computed by integration along the balance line of the zone-1 flux per 
unit length from s = U to s = T. This flux may be expressed as 


P (2) - [ 

t' 1 '* Ju 


t-0) 


W _ ,(») 

1 + 1/2 S l - 1/2 

-f«,( 2 ) _.(*) | 

_ l S / + l/2 S j~ 1/2) 


ds 






(2.5) 


J=P 


s (1) -s (1) 1 

’i+1/2 S l- 1/2 


where the iVy,* values are interpolation coefficients that note the dependence of the flux at 
the jth point (zone 2) on the /th point flux (zone l) and 

l=P 


Once the necessary flux from the adjacent zone is determined, the zonal-boundary 
points of zone 2 may then be updated by solution of equation (2.3). Values of the dependent 
variables in zone 1 at the zonal boundary are obtained by linearly interpolating the newly 
obtained zone-2 solution along the boundary, thus ensuring a continuity of the solution 
across the zonal interface. 

Global conservation is maintained identically in the discrete sense with the above 
procedure. Flux conservation across the flux-balance line CD (see fig. 2) requires the 
satisfaction of the following equation (from Rai, 1986) 

JMAX - 1 

5(*f”+*fS„> + E 

3 = 2 

, LMAX - 1 

=^'' 1, + ^E*)+ e V" 

1=2 

(As before, j points are zone-2 points extrapolated to the balance line and l points are zone- 
1 points along this line.) As noted by Rai (1986), the left and right sides of equation (2.6) 
are but discrete forms of the line integral of the numerical flux F along the flux-balance 
line with zone-2 and zone-1 points, respectively. The equating of the two discrete integrals 
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represents conservation of flux across the balance line and hence the zonal boundary. 
Equation (2.6) does not uniquely define the unknown Pj 2 \ however. Equation (2.5) may 
be shown to satisfy the above statement of global conservation while uniquely specifying 
these fluxes. 

The zonal scheme is not restricted to problems with simple, two-zone interfaces as in 
figure 2, but instead may be used on patched grids with an arbitrary number and placement 
of zonal boundaries. Consider, for example, the six-zone grid for the computation of flow 
about a double-airfoil configuration (fig. 4). This grid is complicated by the intersection 
of zonal boundaries (e.g., point A). Zonal-interface points near these intersections require 
careful attention to strictly maintain conservation; often the flux cells associated with these 
points have six or more sides. 

Treatment of these unusual cells may be explained with reference to figure 5. The 
flux cells depicted here (denoted by dashed lines and bounded on at least one edge by a 
flux-balance line) are typical of those found near the zonal-boundary intersections in the 
grid of figure 4, where zones 1, 2, and 3 or zones 4, 5, and 6 meet. Note that for the 
given choice of flux-balance lines, one zonal-boundary cell is constrained to be six-sided. 
(Points along zonal boundary AA in zone 3 and zonal boundary BB in zone 1 are updated 
by the integration of eq. (2.3) using the flux cells as shown, whereas points along these 
boundaries in zone 2 are interpolated from this solution.) To update point K of zone 3, 
the £y +1 / 2) A, Ej- i/ 2 ,fc> and Fj t k-i /2 fluxes are obtained as usual, using the cell sides EF, 
BA, and FA, respectively, with metrics such as 

Zt) = ( x)e - ( z)f for Ej + i/2,k 

Xi = ( x)f - {x)a for Fj,k -\/2 


* A 

However, must be modified to E*j_ l ^ 2 k to include the flux through side DC by 

adding an integration along the flux-balance line from C to D. 


^-w-A-w + / c 


Ei 


S/- 1 / 2 ] 


ds 
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Similarly, the flux Fj^+ 1/2 is the sum of two terms, each term an integration of fluxes 
along a different flux-balance line 


Fj,k+ 1/2 - J 


: ds 


+ 


B [ 5 m+l/2 ~ s m-l/2 

* F n 


J D Z 5 n 


D [^n+l/2 *®n— 1 / 2 ] 


ds 


(The subscripts /,m, and n are running indices along the flux-balance lines.) Given these 
flux definitions, point K may now be explicitly updated using equation (2.3). 

The patched grid approach allows independent grid generation in each zone. This 
flexibility may result in meshes with discontinuous grid lines or grid-line slopes or both at 
zonal boundaries. Slope or metric discontinuities require attention for accuracy reasons. 
Consider the points along the symmetry line in the grid of figure 4. One such flux cell 
is sketched in figure 6. Although there is grid-point continuity between the two zones 
(and therefore no need for special flux interpolations), the use of the standard central- 
differencing technique for evaluating the metrics in the ^-direction leads to inaccuracies. 
(Depending on the chosen numerical method, numerical flux may be a nonlinear function 
of the slope of the grid lines.) A recommended treatment of these cells is to consider 
them as six-sided and to compute fluxes through each face, using metrics appropriate to 
each side (BA, CB, DC, DE, EF, and FA). This idea is also used in conjunction with 
point-discontinuous boundaries. 


2.3 Summary of Scheme Properties 

The results presented by Rai (1986) and in other publications (Hessenius and Rai, 
1986; Rai, 1985a; Rai, 1985b; Rai, 1985c) illustrate the favorable properties of this patched 
grid scheme. Solution discontinuities, such as shock waves and contact surfaces, will move 
smoothly through grid interfaces with minimal distortion to assume their proper location 
and strength, a direct consequence of the conservative nature of the scheme. Solution 
contours are smooth and continuous across the zonal boundaries. The procedure is nu- 
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merically stable despite the transition of strong discontinuities from zone to zone. Also, 
vortices are shown to move freely through zonal interfaces without distortion (Rai, 1985). 

As reported by Rai (1986), the zonal-boundary flux interpolation (eq. (2.5)) is not 
perfectly free-stream preserving. However, the drift in the dependent variables near the 
zonal boundary when integrating free-stream conditions (without boundary conditions) is 
due to a term that is second-order in magnitude and is proportional to the curvature of 
this boundary. Hence, the departure from free stream may be lessened by increasing the 
number of grid points along a curved zonal boundary. Rai reports free-stream drifts in 
density of only 0.1% after repeated integrations on a grid such as that of figure 2. 

The references cited in the preceding paragraphs describe extensions of this two- 
dimensional method for use with second-order accurate, implicit integration schemes. The 
zonal boundary condition can also be made time-accurate, a necessity for predicting flow 
fields in which zones move with respect to one another such as in a zonal computation of 
flow about a helicopter rotor/fuselage. 
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2 Two-zone grid illustrating zonal scheme in curvilinear coordinates 





Six-zone grid for a double airfoil configuration 
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CHAPTER 3 


APPLICATIONS OF RAI’S TWO-DIMENSIONAL 
PATCHED GRID SCHEME 


Rai’s results (1986) clearly indicate the feasibility of using a conservative grid-patching 
scheme in the numerical simulation of fluid flow. This chapter contains a validation of Rai’s 
zonal blunt body computation and a first demonstration of the usefulness of this patched 
grid method in the solution of geometrically complex problems. 


3.1 Zonal Blunt-Body Computation 

To substantiate the results of Rai (1986) and to validate a newly developed two- 
dimensional Euler code, the supersonic flow about a blunt body (as depicted in fig. 2) 
was computed at a free-stream Mach number of 2.0. The numerical procedure used at 
interior points and at the zonal boundary was described in Chapter 2. Equation (2.3) 
was integrated to convergence using the first-order accurate, explicit, split-flux scheme of 
Steger and Warming (1981). The numerical fluxes E and F for the Steger and Warming 
scheme may be written as 


Ej+ 1/2,* = E-(Q i+1 , k ) + E + (Q,, t ) 

Ej,k+ 1/2 = E (Qy,k+i ) + F + (Qj,k) 


(3.1) 


where 



and Q, E, and F are defined in equation (2.2). The flow field was initialized with free- 
stream conditions everywhere. Flow tangency was enforced at the cylinder boundary, and 
a reflection principle was used to enforce symmetry along the lower boundary (see fig. 
2). The free-stream boundary was held fixed throughout the computation; the supersonic 
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exit boundary was naturally treated with purely backward ^-differences by the split-flux 
formulation. 

Computations were performed on the patched, two-zone grids of figures 7a and 8a. 
The zonal interface in the grid of figure 8a was designed to coincide with the converged 
shock location. In the transient solution, a bow shock forms at the body and moves leftward 
(through the zonal interface in fig. 7a) to its steady-state position. The converged results 
of figures 7b and 8b (pressure contours) are typical of first-order accurate computations of 
this flow. As in Rai’s calculations for this same flow, the captured shock lies slightly to the 
left of the converged shock location of Lyubimov and Rusanov (1973) (shown as square 
symbols in these figures). Note the smoothness of the contours across the zonal boundary 
despite the discontinuity in grid lines. No convergence problems were found with either 
computation. These results are of identical quality to those of Rai (1986). 

3.2 Applications to Geometrically Complex Problems 

The ability of the patched grid method to simplify grid generation about complex 
topologies, and in doing so, to easily enable grid refinement in selected regions is demon- 
strated in the zonal solution of a multiple-body problem. A four-zone solution was com- 
puted about a double-airfoil configuration (both NACA 0006 airfoils); the grid is shown 
in figure 9. With the addition of appropriate source terms in the governing equations for 
axisymmetric flow, this calculation could represent flow about a hollow projectile or an ax- 
isymmetric cowl, as is the case studied by Nietubicz et al. (1979). (The results presented 
here are two dimensional, however.) Unlike that of Nietubicz et al., grid generation for 
the problem was a relatively simple task. An elliptic grid generator was used to obtain the 
meshes in zone 1 and zone 3 (zones 2 and 4 are merely reflections of 1 and 3, respectively). 
The use of a C mesh about the airfoils avoids the introduction of an H-grid singularity 
at the leading edge. The subregional meshes were then patched together to form a global 
grid. Note the discontinuity of grid lines along the zonal boundaries between regions 1 and 
3 and regions 2 and 4 (fig. 9b). 
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The computational results for the flow about the two-airfoil configuration were ob- 
tained by integrating equation (2.3) to convergence using the first-order-accurate, explicit 
Osher scheme with numerical flux defined as (from Chakravarthy and Osher, 1982) 



As in the blunt-body computations, the dependent variables were initially set equal to free- 
stream values. The inflow and upper/lower boundaries were fixed at free-stream conditions. 
Tangency was enforced along the airfoil surfaces, and the supersonic flow at the downstream 
exit boundary was properly treated with backward differences of the governing equations. 

Converged pressure contours are presented in figure 10 for the flow about the double- 
airfoil configuration at free-stream Mach number of 1.5 and zero angle of attack. The two 
bow shocks upstream of the airfoils are apparent in figure 10. These shocks intersect in 
the region between the airfoils and are very quickly dissipated thereafter by the expansion 
waves that are generated by the surface curvature of the airfoils. The flow then expands 
to a pressure minimum in this “nozzle” region and approaches free-stream conditions at 
the outflow boundary. Note that the contours are continuous across the zonal boundaries 
(shown as dashed lines), even in shock regions. The smoothness of the contour lines is 
a result of the conservative nature of the scheme and the enforcement of continuity in 
interpolating values at the zonal boundary. 

The choice of the grid system shown in figure 9 is an arbitrary one. It is possible to 
zone a flow field in a number of different ways and obtain results, but as in any compu- 
tation, there are usually advantages of one grid over another. The grid system of figure 9 
is preferred for this computation. The C mesh provides the best treatment of the airfoil 
trailing-edge, the wake, and the nozzle region because of the orthogonality of the grid. This 
system has additional advantages in that it minimizes the number of zones and it provides 
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simple, two-zone interfaces. Nevertheless, one may wish to use a different patched-grid 
system for the same problem, such as that of figure 11. Again, the airfoils are contained 
within body-conforming meshes (this time O grids) which are easily generated with an el- 
liptic grid generator (GRAPE, from Sorenson, 1980). The meshes in the remaining regions 
are simply generated algebraically. Pressure contours for this grid are shown in figure 12 
(identical contour levels are plotted in figs. 10 and 12). The same flow-field features are 
apparent in figure 12, with a pressure maximum behind the bow-shock intersection and 
an overexpansion in the nozzle region before returning to free-stream conditions. Again, 
the contour lines pass smoothly through the zonal boundaries. (An occasional mismatch 
of contour lines may be observed at zonal boundaries. The interpolation technique in the 
contour-plotting routine is different from the interpolation procedure used at the zonal 
boundary; hence, small errors may arise in plotting.) 

As expected, only qualitative comparisons may be made between the first-order results 
of the two zonal grid systems. Close quantitative agreement would require second-order- 
accurate computations which are less grid dependent. The six-zone results (fig. 12), 
although believed to be inferior to those shown in figure 10 because of the metric disconti- 
nuities and cell skewness in the grid of figure 11, are included to demonstrate the capability 
of the zonal scheme to produce stable and smooth solutions given an arbitrary zoning of 
the flow field. 

3.3 Other Possible Two-Dimensional Applications 

The difficulty in generating a single grid for regions about complicated configurations 
is alleviated with the patched grid technique. The use of this technique is now described 
in the case of two other complex, two-dimensional problems of current interest. 

The flow about a tri-element airfoil, or “augmentor wing” airfoil, has been computed 
by several researchers (Lasinski et al., 1982; Flores et al., 1984). Although a type of zoning 
was employed by Lasinski et al. (1982), the coordinate lines of the grid were still required 
to be continuous across zonal interfaces. This restriction resulted in the introduction of 
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grid singularities in the flow field (also the case in Flores et al., 1984), as well as a reduction 
in clustering and spacing control. The ability of the patched grid method to interface grids 
that are discontinuous across zonal boundaries permits the use of a grid system such as 
that of figure 13. The grids in each zone can be individually generated using algebraic or 
elliptic grid generators. Zones may also be refined at will, with no continuity restrictions 
at zonal boundaries. 

Another candidate for the zonal approach using discontinuous grids is shown in figure 
14. The dual-throat nozzle problem (further described by Chakravarthy, 1981) requires a 
fine mesh in zone 2 to adequately resolve a normal shock. The zoning technique allows 
for an abrupt refinement in this region. Note that the zonal grid for this problem may be 
obtained using a simple algebraic procedure. 
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(b) 

Fig. 7 Results of a two-zone, blunt-body computation 

(a) two-zone grid (b) pressure contours at = 2.0 
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Fig. 9 Four-zone grid for the two-airfoil configuration 
(a) zones 1-4. 
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Fig. 9 Four-zone grid for the two-airfoil configuration 
(b) enlarged view of zones 3 and 4 
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Fig. 10 Pressure contours for supersonic flow about the double-airfoil 

configuration, computed on the four-zone grid = 1.5, a = 0° 
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Fig. 11 Six-zone grid for the two-airfoil configuration (a) zones 1-6 
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Fig. 11 Six-zone grid for the two-airfoil configuration 
(b) enlarged view of zones 3 and 6 
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Fig. 13 Zoning suggestion for the “augmentor wing” airfoil configuration 
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Fig. 14 Zoning suggestion for the dual-throat nozzle problem 
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CHAPTER 4 


DEVELOPMENT OF A THREE-DIMENSIONAL 
PATCHED GRID SCHEME 

As researchers in the field of computational fluid dynamics (CFD) attack problems 
of increasing geometric complexity, especially in three dimensions, the need arises for 
a more versatile approach to grid generation than the traditional, single-module type of 
mesh definition. Consider, for example, the discretization of the flow field about the Space- 
Shuttle orbiter with external tank and rocket boosters (fig. 15a). Multiple-body problems, 
such as this one, are challenges using even the most advanced grid-generation techniques 
given the necessary constraints on grid quality (e.g., smoothness, desired clustering in 
gradient regions). 

The zonal concept simplifies mesh generation by partitioning the flow field into sev- 
eral regions within which grids are independently generated. Consider again the problem 
depicted in figure 15a and imagine a partitioning of the flow field at cross-section A as in 
figure 15b. A division of the field into well-shaped, four-sided zones (actually, six-sided 
blocks for this three-dimensional problem) permits the use of existing grid generation tech- 
niques within each zone. The subregional grids could be simply patched together along 
the specified interfaces. A proposed zoning for cross-section B is shown in figure 15c (this 
zoning is completely independent of the one shown in fig. 15b, however). 

Chapter 2 is a review of a conservative, zonal scheme for patched grids in two dimen- 
sions, that is, grids that interface along a curve. (The grid lines may be point-discontinuous 
at the zonal interface, as shown previously in fig. lc.) The procedure has been tested in 
a variety of settings, including the computation of blast- wave diffraction, flow about a 
double-airfoil configuration, and flow about rotor-stator combinations at both subsonic 
and supersonic free-stream Mach numbers. Calculations have demonstrated the ability of 
shocks as well as vortices to pass through zonal interfaces with minimal distortion (Rai, 
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1986; Hessenius and Rai, 1986; Rai, 1985a; Rai, 1985b; Rai, 1985c). 

This two-dimensional patching technique may be readily applied to three-dimensional, 
parabolic problems where grid generation is performed in two dimensional planes (normal 
to the body axis), and the solution is obtained by marching downstream in the streamwise 
coordinate. But a fully three-dimensional, zonal computation using patched; grids requires 
a a division of the flow field into well-proportioned, six-sided blocks. There are three 
levels of difficulty associated with the implementation of a conservative, three-dimensional 
scheme for patched grids. 

1) The simplest form of patching in three dimensions occurs when a coarse grid is refined 
in certain regions by the subdivision of coarse grid cells into several fine grid cells. 
The flux-balancing procedure to ensure conservation at the zonal interface is greatly 
simplified when integer mesh-cell ratios (zone to zone at the interface) are present. 

2) A more general form of patching occurs when arbitrary point distributions are per- 
mitted at the interface in each zone. At this level, however, the interface is restricted 
to be planar in physical space. 

3) The most general form of patching is interfacing along an arbitrary (nonplanar) surface 
with arbitrary point distributions at the zonal interface in each zone. 

The focus of this work is on the development of a three-dimensional, globally conserva- 
tive, boundary scheme for patched grids, applicable in generalized coordinates for arbitrary 
point distributions on a planar zonal surface (level 2 above). As was done in the original 
two-dimensional work by Rai (1986), the three-dimensional scheme is derived within the 
framework of an explicit, first-order method for the solution of the Euler equations. It 
can then be generalized for second-order accurate, implicit computations as it was in two 
dimensions (Rai, 1985). 
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4.1 Three-Dimensional Governing Equations 

The unsteady Euler equations in three dimensions may be written as 


Qt + E x + Fy + G z = 0 (4-1) 

where 
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V = (l ~ l)[e “ \p { u2 + v 2 + ty2 )l 

Again, p is the density; p is the pressure; u, v and w are the velocity components in the x, 
y, and z directions, respectively; and e is the total energy per unit volume. Applying the 
independent- variable transformation 


t =t 


Z = Z{x,y,z,t) 

ri =n(x,y,z,t) 

$ = ((x,y,z,t) 

to equation (4.1), one obtains the Euler equations in generalized coordinates 


Qt + + Fff + Gj = 0 (4.2) 

where 

Q = Q/J 

E(Q , Z) = (ZtQ + ZxE + £ y F + ZzG)/J 
HQ . v) = in tQ + *izE + ri y F + ^G)/ J 

G(Q, f) = (ftQ + t xE + $ y F + $ Z G)/ J 
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with the following metric identities 


J 1 = x i y„z i + x^y^Zr, + x r) y^z i - x^z,, - x„y ? z ? - x i y r ,z i 
lx = J{yr,z< - z„y f ) r\ x = J(y f z ? - z t y ( ) 


( y ~ x n 2 f) 

£z = J{ x t iy$ ~ y»jXf ) 

?x = J{yiZr, - ny r>) 

fy = J ( z t x n - x i z r,) 


ri y = J{z i x i -x^z i ) 
r) z = J(x ? y^ - y ? x $ ) 
ft — ~X t £x ~ yr£y ~ Zrtz 
m = -X r T/i - y T »?y - ZtVz 


fz — J(x^yr) — y$ X Tj) ft — X T f z yrfy Xrfz 

An explicit, conservative, finite-difference scheme for integrating equation (4.2) is given by 


Q n+ ' - q" - ■ E r.-i/3,>,n 

Ar Af 

rn _ E>n An _ An 

J (ij + l/2,*) ^(ij-l/2,fc) _ + 1/2) 1/2) _ 

Ay Af 


(4.3) 


where -E, n +1 / 2 > k > ^tj+i/2,fc an ^ ^T,y,fc+i/2 are numerical fluxes consistent with the trans- 
formed fluxes E , F, and G, respectively. The expressions for the numerical fluxes depend 
on the chosen numerical integration method (discussed in section 4.5). 


4.2 Metric Computations for Free-Stream Preservation 

The computation of numerical flux requires an evaluation of the following metric 
quantities: fi/J, (y/J, Zz/J for E (see eq. (4.2) definition). A proper determination of 
the metrics is necessary for the maintenance of free-stream conditions in undisturbed areas 
of the flow. In free-stream regions, the governing equations degenerate to derivatives of 
metric quantities alone. For example, the continuity, x-momentum, and energy equation 
become, at free stream, the single expression 

Mj) + a,(j) + 3dj) = ° 

(Similar equations containing metric derivatives only can be written for the y- and z- 
momentum equations at free stream.) In three dimensions, simple three-point central 
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differences in the calculation of the inverse metrics (i.e., y ^ , etc.) to construct the 
needed metrics will not guarantee free-stream preservation on an arbitrary grid; that is, 
such a metric construction will not satisfy the above equation within acceptable numerical 
error bounds. The problem arises from the fact that analytical derivative operators are 
commutative, whereas numerical difference operators are not. In two dimensions, it is 
fortuitous that central differencing of the inverse metrics will yield perfect preservation of 
free-stream conditions. 

Pulliam and Steger (1978) suggest two possible approaches to metric computation in 
three dimensions to gain free-stream preservation: 

1) employ special weighted averages in the differencing of the inverse metrics such that 
perfect error cancellation occurs when metrics are differenced in the governing equa- 
tions at free-stream conditions 

2) reformulate the governing equations by subtracting the free-stream fluxes from the 
original equations. 

Kordulla and Vinokur (1983), however, recommend an approach to the problem of met- 
ric definition that is consistent with our “control volume” method of solution at zonal- 
boundary points. In fact, their technique for defining the transformation metrics is stan- 
dard in finite-volume computations. Their procedure is described in the following para- 
graph and is adopted in this work. 

Consider the computational cell about the point (ij,k) in figure 16 (cell sides are 
located halfway between neighboring grid points). The £?(,+i/ 2 ,y,fc) flux is the flux through 

A 

the nonplanar cell face, ABCD. The computation of E requires the metrics £ x / J, and 

AAA 

iz/J which can be interpreted as the *, j, and k components, respectively, of the surface 
vector S A bcd ■ (The unit vectors t, j, and k are along the x, y, and z axes, respectively.) 
The surface vector, 5, is a vector unique to a surface described by four vertices. The 
magnitude of S is an approximation of the area of the surface, and the direction of this 
vector is normal to both diagonals of the surface. The derivation of various forms of 
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the expression for S is given by Kordulla and Vinokur (1983), but the simplest form is 
presented here. 

Sabcd = 2 (^> - *b) x (rc ~ r A ) (4.4) 

In the above expression, ?k is the position vector to point K from an arbitrarily placed 
origin. The components of S a bcd are the metrics £ Z /J, £ y /J, and £ Z /J needed to define 
E through cell-face ABCD. The use of equation (4.4) to define the metrics at all cell 
faces will ensure the preservation of free-stream quantities when integrating the governing 
equations without boundary conditions given free-stream initial conditions. This claim has 
been verified numerically in the computations presented in this work. 

An exact evaluation of the cell volume, J _1 , needed to advance the solution in time, 
is not necessary in steady-state applications since it, in fact, functions only as a scaling of 
the time step. For the steady-state results of Chapter 5, J -1 is simply evaluated at the 
cell center (grid point) using central differencings in the expression of equation (4.2). A 
more accurate computation of the cell volume is necessary for time-dependent calculations; 
interested readers are referred to the work of Kordulla and Vinokur (1983) for one possible 
approach. 

4.3 Solution at Interior Points 

Consider a two-zone, three-dimensional computation such as that depicted in figure 
17a (problem is shown as a quadrilateral in computational (£,»?,?) space, with maximum 
grid dimensions IM^ m \ J and KM^ m ^ where m corresponds to zone number). The 
zonal interface is shown as a shaded plane at tj — JM ^ in zone 1 and rj = 1 in zone 2. 
The interior points (away from boundaries, including zonal boundaries) in each zone are 
explicitly updated by solving equation (4.3) using the procedure in section 4.2 to define 
the metrics. 
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4.4 Solution at the Zonal- Boundary Plane 

4.4.1 Development of Conservative Flux-Interpolation Scheme 

In this section, the development of a conservative flux-interpolation scheme for three- 
dimensional patching at a planar zonal boundary is outlined. The specific details of its 
implementation and coding are discussed in a subsequent section (4.4.2). 

The solution procedure at the zonal boundary is a two-step process. 

1) The zonal-boundary points in one of the two zones are updated using a scheme that 
guarantees global conservation across the boundary. This scheme requires integrating 
the governing equations using equation (4.3) while interpolating for fluxes needed from 
the adjacent zone. 

2) The solution at points on the zonal boundary in the neighboring zone is found by 
interpolating the solution found in step 1 along the zonal interface thus ensuring 
continuity of the dependent variables. 

Consider again the problem shown in figure 17a. A view of the zonal boundary plane 
with f f = constant lines is shown in figure 17b where (x) points belong to zone 2 and 
(•) points belong to zone 1. If one chose to integrate the boundary points in zone 2 while 
interpolating for the solution in zone 1, then the flux cell associated with the point (i2,l,k2) 
on the zonal boundary in zone 2 is like that of figure 18. Notice that the cell is formed 
by extrapolation of grid lines into zone 1 to the flux-balance plane (q = 1/2 in zone 2 or 
V = J — 1/2 in zone 1). Five of the six fluxes necessary to update point (i2,l,k2) 

using equation (4.3), namely, ^ 2 + 1 / 2 , 1 , * 2 )’ ^( 12 - 1 / 2 , l,*2)’ ^(W,3/2,fc2)’ ^(* 2 , 1 , * 2 + 1 / 2 )’ and 
^(i 2 ,i,Jfc 2 — 1 / 2 )’ are readil y formed from values of the dependent variables in zone 2. The 
flux through the flux-balance plane, -^ 21 / 2 * 2 )’ * s obtained by interpolation of zone-1 
fluxes, as explained in the following paragraphs. 

A planform view of the flux-balance plane is presented in figure 19 showing the pro- 
jections of both zone 1 and zone 2 grid points on this plane. The locations of zone 1 points 
(•) on the plane are found by simply averaging the coordinate values at r/ = JM ^ and 
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t) = JM — 1. The locations of zone 2 points (x) require extrapolation of zone 2 grid lines 
to this imaginary j2=l/2 surface (a detailed description of the extrapolation is given in 
the next section). The cell face associated with the point (i2,l,k2) and the flux 1/2 Jt2) 
is labeled as RSTU. Cell faces of zone 1 points are also shown by segmented lines. 

The zonal-boundary scheme developed for two-dimensional problems in Rai’s work 
(1986) and for three-dimensional applications as described in this section guarantees global 
conservation of mass, momentum, and energy in a discrete sense at the zonal interface. To 
affect this global conservation, the surface integral of the zone 1 numerical flux, F, on the 
flux-balance plane must be equal to the surface integral of the zone 2 numerical flux on 
this same plane. The discretized form of this requirement may be written 
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where jl = - 1/2 and j2 = 1/2. The left and right sides of this equation are the 

discrete forms of the surface integral of the numerical flux, F, on the flux-balance plane in 
zone 1 and zone 2, respectively. Like its two-dimensional counterpart (eq. (2.6)), equation 
(4.5) is a necessary, but not a sufficient condition to uniquely define the unknown F^ 2]2 k2 ^ 
fluxes. However, a generalization of equation (2.5) for three-dimensional application does 
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satisfy the above equation of global conservation. Again referring to figure 19 and the flux 
cell face RSTU on this planar surface. 


j>(2) 

r (x2,j2,k2) 


f'f 


T irU) 


I (jU) 


dx dz 


fci) I 


(4.6) 


In essence, the integrand of equation (4.6) is a representation of the zone-1 flux per unit 
area at the flux-balance surface. (A piecewise-constant distribution of numerical flux over 
a cell face is assumed in this work.) Zone-1 flux per unit area is then integrated over the 
area of the zone-2 cell face to yield the total flux through face RSTU. In discrete form, 
equation (4.6) becomes 


KM 11 * IM {1) M 1 ) 

p( 2) _ | c( 2 ) I V' V' Ar(«2,j2,fc2) r (iljl,kl) 

r (x2,j2,k2) ~ ■ (t’2 j’2,fc2) • ZL iV (tl,jl,fcl) , -(1) 

*1 = 1 *1 = 1 P(*l,.jl,A:l) 


(4.7) 


(t2 ?2 fc2) • * (2) 

The ji ’*•!) are interpolation coefficients that represent the dependence of j2 k2 } 

on the j fluxes. Physically, these terms are nothing more than ratios of the areas 

of overlap of RSTU with zone-1 cell faces to the total area of RSTU, that is 


Ai2,j2,k2) 

Ar(*2,j2,fc2) 

iV (*l,yi.fcl) ~ | q(2) 

\°(x2,j2,k2) 


(4.8) 


where denotes the area of overlap of the cell face associated with (il jl,kl) with 

that of (i2j2,k2) on the flux-balance plane. One final condition is necessary to ensure that 
equation (4.7) satisfies the statement of global conservation (eq. (4.5)), that is 


KM il) IM il) 

E V J ^(*2, >2,1:2) 

2L, JV (tl,.)l,kl) 
Jfcl = l *1=1 


= 1 


(4.9) 


Equation (4.9) states that the area of a cell face in zone 2 on the flux-balance plane is 
identically equal to a sum of the areas of zone-1 cell faces, ensuring that the flux from the 
zone-1 side of the flux-balance plane is completely transferred to the zone-2 side of this 
plane. 
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Given the definition of flux from a neighboring zone (eq. (4.7)), it is now possible 
to update all points on the zonal-boundary plane in zone 2 by solution of equation (4.3). 
The solution at points on this plane in zone 1 is then found by interpolating the zone-2 
dependent variables on the interface using a procedure described in Appendix A. 

4.4.2 Implementation of Zonal-Boundary Scheme 

Equation (4.7) defines the flux through a flux-balance plane and therefore the flux 
from a neighboring zone at the zonal boundary. Use of equation (4.7) will, by definition, 
ensure global conservation in the discrete sense. The F^j l kl ) terms in this expression 
are formed in a conventional manner using the definition of numerical flux for a given 
numerical integration scheme. The metric terms (components of S) are computed using 
equation (4.4). The interpolation coefficients, of equation (4.8), are the remain- 

ing quantities to be determined. A schematic of the procedure for computing these terms 
is shown in figure 20 and is discussed in the following paragraphs. 

Assuming the integration of zone-2 points on the zonal-boundary plane (followed by 
an interpolation for values at zone-1 points on this boundary), the flux-balance plane is 
located in zone 1 halfway between the last two constant ^-coordinate planes (see again 
fig. 17a). It is necessary to find the projections of both zone-1 and zone-2 points on 
this ficticious plane. The locations of zone-1 points are found by simply averaging the 
coordinate values at r) = J and r) = — 1. The locations of zone 2 points require 

an extrapolation of zone-2 grid lines to this imaginary j2=l/2 surface. The extrapolation 
is performed in such a way as to emulate the point relationships on the zonal-boundary 
surface using linear interpolation over triangles as described by Chung (1978) (see also 
Appendix A). 

Since only three points define a plane, it is necessary to work with triangular portions 
of each cell face to determine areas of overlap. The cell area on the flux-balance plane 
associated with each zone-2 point is now considered in turn and is divided into a number 
of triangles. On a rectangular, planar zonal boundary, each cell is adequately described 
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using only two triangles. Similarly, the cell areas associated with zone-1 points on the flux- 
balance plane must also be “triangulated.” Given the triangles associated with a zone-2 
point, one must now compute the areas of overlap with all zone-1 point triangles. The 
total area of overlap of all triangles of zone-2 point (i2J2,k2) with all triangles of zone-1 
point (il jl,kl) is stored in the coefficient -A()i 

The problem in computing the flux-interpolation coefficients is now reduced 

to finding the common area of two overlapping triangles. An unsuccessful attempt was 
made to develop an analytical expression for determining the polygon of overlap and hence 
the area. If such an expression is indeed derivable, computation of the jfeij* would 

be inexpensive enough for unsteady (moving grid) applications where a new evaluation is 
necessary at each time step. In this work, however, an algorithmic approach was employed 
that, although computationally intensive, resulted in an insignificant amount of computing 
overhead because only one evaluation of the is needed for these applications. 

A modified Sutherland-Hodgman clipping algorithm (as outlined by Newman and Sproull, 
1979, see Appendix B), borrowed from the field of computer graphics, is used to find the 
vertices of the polygon of overlap (fig. 21). The area of any convex polygon with n ordered 
vertices may be found by the following expression 

1 n 

area = -| ^(xi+iy, - x,!/»+i)| (4.10) 

t=i 

where x n+i = x x and y n +i = yi- (The polygon resulting from the overlap of two triangles 
will, conveniently, always be convex.) The accuracy with which the overlapping areas are 
computed determines the degree of free-stream preservation at a planar zonal boundary. 
In this case, triangular overlaps can be computed to machine accuracy enabling perfect 
free-stream preservation. 

The procedure depicted in figure 20 (and described above) is by no means optimal. 
There are “smart” ways to avoid unnecessary computation by eliminating the triangulation 
of zone-1 points far from the zone-2 point under consideration. 
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4.5 Numerical Integration Scheme: Osher Algorithm 

In the development of the three-dimensional patched grid method in the previous 
section, it is assumed that a conservative finite-difference scheme is used to define the 
numerical fluxes E, F, and G. The results which follow in Chapter 5 were obtained with 
the first-order accurate Osher algorithm. This scheme has been used extensively in two 
dimensions and is outlined for three-dimensional calculations in Osher and Chakravarthy’s 
work (1983). The algorithm as presented in the reference is unsuitable for general com- 
putations, however, since it requires that £ z , r) x , and $ x are nowhere zero. A necessary 
modification of the scheme for an arbitrary grid is described in this section. 

The Osher scheme is an upwind, shock-capturing algorithm based on a Riemann 
solver. Unlike an exact Riemann solver, it uses compression waves to approximate shocks, 
thereby simplifying the algorithm. The numerical flux functions may be written as 


•Qi+l.j.k ( Qg 


rQi.j+i.k / d p\ 

Fij+i/2,k = Fi,j+i,k - l an I 

rQi.j.k + i f dG\ + 

Gi,j,k+ 1/2 = G i,j,k + 1 - ^7) I 

' Qi.j.k \ a{ * ) 


where E, F, and G are given in equation (4.2). The integrals are evaluated along subpaths 
in state space chosen to coincide with the right eigenvectors of the Jacobian matrices 
dE/dQ , dF/dQ , and dG/dQ for the £, rj, and f directions, respectively. Figure 22 
depicts the three subpaths between any two adjacent grid points (m-1 and m or m and 
m+l) where o denotes the eigenvalue associated with each subpath. (Subpath end points 
or states are labeled at 1/3 intervals.) The “-j-”in equation (4.11) indicates contributions 
resulting from positive eigenvalues along the subpaths, i.e., the subintegrals will be zero if 
the corresponding eigenvalue is negative all along the subpath. If not zero, the subintegrals 
are simply flux differences between the end points of the subpath, assuming there are no 
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sonic points along the path. The treatment of sonic points is outlined fully by Osher and 
Chakravarthy (1983). 

Key to the development of the Osher algorithm is the identification of invariant quan- 
tities, called Riemann invariants, along each subpath. By definition, Riemann invariants, 
'k, associated with the nth eigenvalue must satisfy 


Vq* r n (Q)=0 (4.12) 

where r n denotes the right eigenvector associated with the nth eigenvalue. It may be 
shown that, given the above definition, the quantities are invariant along their distinctive 
subpath (Chakravarthy and Osher, 1982). (Equation (4.12) is the only condition necessary 
to qualify a Riemann invariant.) Consequently they will be used to determine values of 
the dependent variables, Q, at the intermediate states and, ultimately, the flux at these 
points. 

Figure 22 presents Riemann invariants suggested by Osher and Chakravarthy (1983) 
for three-dimensional computations. Equating these invariants along their respective sub- 
paths with known values at grid points (for example, m-1 and m) produces ten equations 
for the ten unknown quantities at intermediate states (m-2/3 and m-1/3). Note that 
should £ x , rix, or £ x be anywhere zero, v and u> are no longer linearly independent, and the 
equation set degenerates to nine linearly independent equations for ten unknowns. Hence, 
it is necessary to identify a new v and tv that are always independent along the 1, 4, 3, 
and 6 subpaths. 

Assume the following form for v and w, where c, are constants yet to be determined 

V = Ci« + C 2V + C3W 

(4.13) 

W = C4U -|- C 5 t> + CqW 

Taking the gradient of both v and w with respect to the dependent variables Q produces 
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the row vectors 


Vgt) = 
VqU) = 


,c 1 C 2 C3 x Cl C2 C3 _ 

--«+ —« + -to), -i, -, — , 0 

PPP ppp 

,C\ C5 Ce . C4 C5 Ce 

U + — V + — u>), — , — , — , 0 

ppp ppp 


(4.14) 


The eigenvector associated with the o = D(u — c) subpath (fig. 22) provides the most 
stringent test for the qualification of a Riemann invariant by equation (4.12) and may be 
written as 


k z p--u, k y p - -v, k z p--w , £ 

L c c c c J 


(4.15) 


where 


, = ^ ±^±^l - kzPU - k„pv - k z pw + 


2c ~ (7-I) 

From equation (4.12) (using eqs. (4.14) and (4.15)), the following constraints on the 
constants, c^, are derived 

c\k x + c 2 k y + c 3 k z = 0 

(4.16) 


c 4 k x T c§ky "I - CQk z — 0 

Note that the above constraints are satisfied if the c t are chosen such that the vectors 
[ci,C2,C3] and [04,05, Co] are orthogonal to the particular metric [ k x ,k y ,k z ]. Furthermore, 
the vectors [01,02,03] and [04,05,03] must be linearly independent to ensure the same for v 
and w. One way to meet all requirements is to rotate the [k x , k y , k z ] vector to a coordinate 
axis in xyz-space and choose the other two coordinate axes as the transformed “c”vectors. 
For instance, if R is a rotation matrix and 


then 




c 4 

C5 

«6 


= R' 1 


= R~ l 
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Details of this procedure are given in Appendix C. 

Excluding a redefinition of the invariants v and u>, there are no other required changes 
to the algorithm for use in general computations. The reader is referred to the work of 
Osher and Chakravarthy (1983) for a more detailed description of the scheme, including 
the treatment of sonic points as they occur along subpaths. 


Fig. 15 (a) View of the Space Shuttle orbiter with external tank 

and rocket boosters 
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Fig. 15 (b) Zoning suggestion at cross-section A 
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ZONE 1 
(•) 

(b) 


Fig. 17 Sketch of a two-zone, three-dimensional problem 

(a) computational domain (b) enlarged view of zonal-boundary plane 
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Fig. 18 Computational cell associated with the point (i2,l 5 k2) 
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Fig. 19 Planform view of the flux-balance plane 
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flux-interpolation coefficients, -A(|i 


















Fig. 21 Sketch depicting the overlap of two triangles. 
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m-2/3 a = Du m-1/3 m+1/3 a = Du m+2/3 



R I EM ANN INVARIANTS ALONG SUBPATHS 
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Fig. 22 Schematic representation of Osher algorithm and the 
identification of Riemann invariants along subpaths. 
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CHAPTER 5 


APPLICATION OF A THREE-DIMENSIONAL 
PATCHED GRID SCHEME 

The three-dimensional interfacing technique developed in Chapter 4 for patched grids 
is used in the computation of supersonic flow about a canard/ wing combination mounted 
on a wall. This problem was chosen to demonstrate the use of grid patching in the solution 
of flow about a multiple-body configuration that is typically difficult to grid. The results 
of this computation indicate that a stable solution can be obtained with the conservative 
interfacing technique and that smooth solutions are generated despite discontinuities in 
grid lines. 

5.1 Problem Description 

A planform sketch of the problem is shown in figure 23a. The canard section is a 3% 
biconvex airfoil, the wing section a NACA 0006. Both canard and wing are of constant 
chord with aspect ratios of 1.5 and 2.0, respectively, and are swept at an angle of 25°. 
Single-module type grid generation for this configuration would be difficult because of the 
y-coordinate offset of the wing from the canard (see fig. 23b, offset is 10% of canard chord). 
Hence this problem is a good candidate for a patched grid treatment. 

The sharp leading and trailing edges of the canard suggest an H-grid treatment of 
the flow region, while the blunt leading edge of the wing requires a C-grid. In a two- 
zone, patched grid approach to the configuration, a 49 x 50 x 26 point H-grid about the 
canard (grid dimensions are in the streamwise, nearly normal and spanwise directions, 
respectively) interfaces with a 111 x 30 x 34 point C-grid about the wing along a planar- 
zonal boundary (also swept at 25°). The wall-plane cross section of the grid (an x-y plane) 
is shown in figure 24. Note that the grid lines are point discontinuous at the zonal interface 
in the y-direction. (The wing grid has been clustered in the streamwise direction near the 
zonal interface to match the streamwise spacing in the canard grid, thus avoiding metric 
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discontinuities arising from abrupt changes in spacing.) Grid-line discontinuities also exist 
in the spanwise or z-direction on the zonal-boundary plane. In fact, there are only three 
spanwise planes that are coincident in the canard and wing zones: the wall plane (minimum 
z), the domain boundary plane (maximum z), and a plane through the canard tip. The 
constant-section, canard and wing grids transition at the tip of each body to grids about 
flat plates with zero thickness. A view of the grid along the span of the canard (a y-z 
plane) is shown in figure 25. (The wing grid in this region is identical in nature.) Figure 
26 presents a three-dimensional perspective of the wall-plane grid with the canard and 
wing. 

The canard and wing zones are shown as quadrilaterals in computational space in 
figure 27. The boundary conditions are indicated on the faces of each computational 
domain. Tangency is enforced on body surfaces. However, there is flow through the 
flat-plate extensions of the wing and canard, as well as through the wake of the wing. 
A symmetry condition is imposed at the wall plane, and conservative flux interpolations 
are performed on the zonal-boundary plane. The inflow is fixed at supersonic free-stream 
conditions. All other boundaries in the flow domain are sufficiently distant to maintain free- 
stream flow with the exception of the downstream boundary. This boundary (supersonic 
outflow in the wing domain) is updated by extrapolating from the interior. 

Although geometrically a two-zone problem, the flow field was actually divided into 
seven blocks (four blocks in the spanwise direction of the wing and three in the spanwise 
direction of the canard, as shown in fig. 28). This blocking was performed to circumvent 
in-core memory limitations of the Cray X-MP. Therefore, only one block of data resides in 
core at a given time while the remaining data resides in easily accessible extended storage 
(solid-state device). The blocks, which are processed sequentially, overlap by one plane in 
the spanwise (f) direction to readily permit the transfer of data between blocks. (Since 
the grid lines are continuous from spanwise block to spanwise block, it is not necessary to 
use the patched grid scheme at these block interfaces.) 
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For this application, the entire flow field is initialized with free-stream conditions at 
a Mach number of 2.0. Equation (4.3), with the first-order Osher scheme flux definitions 
of equation (4.11), is integrated to convergence. The solution is considered converged if 
the maximum change in the density (from step to step) divided by the time step does 
not exceed 10 -5 . This convergence criterion is met after approximately 8000 time steps 
on the grids previously described. The canard-wing calculations presented here required 
approximately 24 hours of Cray X-MP computer time using a fully vectorized version of 
the flow solver. The advantage of extending this zonal-boundary procedure for use with 
implicit schemes (with their characteristic accelerated convergence rates) is evident. 

5.2 Computed Results 

Figures 29 and 30 present the results of a computation of the flow at a = 0°. (In 
all figures the flow is generally from left to right). Figures 29a and 29b depict pressure 
contours in both zones on a plane through the canard tip (one of the three coincident 
z-planes). Note the expected symmetry of the canard contours (from top to bottom), as 
well as that of the attached leading- and trailing-edge shocks (weakened near the tip). The 
trailing-edge shock waves pass through the zonal-boundary plane to intercept the detached 
wing-bow shock. A very slight asymmetry in the wing solution results from the offset of 
the canard chord from the wing in the y-direction. Despite the grid-line discontinuities 
in the y-direction and the z-direction on either side of this spanwise plane, the solution 
near the zonal boundary between canard and wing is smooth and continuous. As in the 
two-dimensional applications, an occasional mismatch of contour lines may be observed at 
the zonal plane. The interpolation technique in the contour-plotting routine is different 
from the interpolation procedure used at the zonal boundary; hence, small errors may arise 
in plotting. 

Pressure contours on the upper surfaces of both canard and wing are displayed in 
figure 30. Again there is evidence of leading- and trailing-edge canard shocks, strong 
gradients near the leading edge of the wing, and a wing trailing-edge shock. The flow over 
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the canard is primarily two dimensional except near the tip within the Mach cones (see 
again fig. 23a, n is the Mach angle). Likewise, the wing experiences little departure from 
two-dimensional flow except at the wall and at the tip. 

Figures 31 and 32 present results for the flow over canard and wing at a = 2°. Again, 
pressure contours are shown in the spanwise plane through the canard tip (fig. 31). At 
a = 2°, there is a leading-edge expansion on the upper surface of the canard with a lower- 
surface shock and a decided asymmetry in the wing solution with greater expansion on 
the wing upper surface as compared to figure 29. The contour lines indicate that the 
solution is smooth and continuous in the vicinity of the zonal boundary. In figure 32, 
note that the wing upper-surface pressures are significantly altered in the spanwise region 
near the canard tip because of canard downwash. Although in theory one would expect to 
observe a relatively weak tip vortex formation at a = 2°, no vortical motion is detected 
in these results. The tracking of such a vortex requires at least a second-order accurate 
numerical integration scheme in conjunction with a very fine grid. The dissipative nature of 
a first-order accurate numerical method precludes the capture and propagation of vortical 
motion. 

As a form of code validation, the computed canard pressures at a — 0° and a = 2° flow 
conditions are compared with linear theory (a numerical solution of the Prandtl-Glauert 
equation from the code PAN AIR as described by Carmichael and Erickson, 1981) in figures 
33 and 34. (Since the flow is supersonic at the zonal-boundary plane and therefore there is 
no upstream influence of the wing on the canard, it is appropriate to compare the present 
results with theory for a canard alone.) The two results are in good agreement throughout 
the span. The slight discrepancies in pressure noted at the leading edge are the result of 
the smearing of the leading-edge shock by the first-order numerical method of the present 
computation. 
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5.3 Discussion 


The canard/wing application was chosen to demonstrate the usefulness of a grid- 
patching technique in the solution of geometrically complex problems. The intention of 
this work is to provide a three-dimensional extension of Rai’s (1986) two-dimensional 
zonal-grid scheme that is globally conservative and stable. The preceding computations 
indicate that this goal has been reached. The calculations are conservative by the very 
nature of the flux interpolation scheme developed in Chapter 4. By definition, the fluxes 
at the zonal boundary satisfy a statement of discrete conservation. Other researchers 
(Benek et al., 1983; Rai, 1985c) working with overlapping grid systems suggest that the 
nonconservative nature of their interface schemes gives rise to high-frequency oscillations 
at the zonal boundary. The solutions presented here (as well as the second-order accurate 
results of Rai (1985) in two dimensions) are free from oscillations in the zonal-boundary 
region. 

The ability to obtain a stable solution using grid patching, however, requires some 
care in grid construction. In performing stability analyses for various two-dimensional 
interface schemes, Eriksson and Rai (to be published) find that if using a central-difference 
integration method, it is necessary to integrate the points on the zonal boundary of the 
coarser grid and interpolate for the points on the boundary of the finer grid. An instability 
will result if coarse grid fluxes are divided to form fluxes for the finer grid using Rai’s (1986) 
approach to grid patching in two dimensions. Although Eriksson and Rai did not observe 
such behavior with upwind integration methods (given modest grid-density ratios across 
the boundary), their work suggests that the choice of zonal-boundary points to integrate 
versus interpolate is not as arbitrary as once thought. 

The present work supports this conclusion. The first grid systems generated for this 
problem were characterized by up to 10 to 1 grid-density ratios from zone to zone in 
the spanwise direction. Points were clustered in the region of the canard tip for body 
definition, but there was no corresponding resolution in the wing zone interfacing with 
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this region. Although the coarser (wing) grid was not expected to support the regional 
gradients as accurately as the finer (canard) grid, no stability or convergence problems 
were anticipated. However, when the coarse grid fluxes were divided to form fluxes for the 
finer grid (as would be the case if the wing zonal- boundary points were integrated), the 
computation would not converge while errors (build-up of pressure contours) accumulated 
at the zonal interface. The stable results presented here were obtained by both reducing 
the grid-density ratio to at most 3 to 1 and also by changing the zone of integration such 
that the flux for a coarse grid cell was specified by the summation of appropriate fine grid 
fluxes. The unreasonable matching of a fine grid with a coarse one in a gradient region 
has an expected effect on solution accuracy but also introduces a stability problem. For 
this reason, patched grid systems should be characterized by modest grid-density ratios at 
the zonal boundary (less than 3 to 1 in gradient regions). Furthermore, the coarser side of 
the interface should be chosen as the zone of integration versus interpolation. 

The numerical integration scheme used in this study, the Osher upwind method, is 
first-order accurate. Likewise, the assumptions made in the construction of the zonal- 
boundary fluxes (piecewise-constant variation of fluxes) and the use of linear interpolation 
to enforce continuity at the patch boundary produce first-order accurate results at the 
interface. 
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Fig. 24 Two-zone grid at wall-plane cross-section (minimum z). 
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Fig. 25 Spanwise view of grid at canard tip 
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Fig. 26 Three-dimensional view of wing-canard and wall-plane grid 
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Fig. 27 Boundary conditions in the computational domain 
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Fig. 28 


Sketch depicting the blocking of a computational domain 



(a) view of entire flow field. 



(b) close-up of zonal boundary. 


Fig. 29 Pressure contours in the plane of the canard tip, M = 2, a = 0° 
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Fig. 30 Upper-surface pressure contours on canard and wing, 

A/oo = 2, a = 0° 


Pressure 



(a) view of entire flow field. 

Fig. 31 Pressure contours in the plane of the canard tip, M 0 Q — 2, a = 2° 
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(b) close-up of zonal boundary. 


Fig. 31 Pressure contours in the plane of the canard tip, M 0 0 = 2, o = 2 C 



Fig. 32 Upper-surface pressure contours on canard and wing, 

Moo =2, a = 2° 
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Fig. 33 Canard pressure coefficient comparisons with theory at 
various span stations, M 0 0 = 2, a = 0° 
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Fig. 34 Canard pressure coefficient comparisons with theory at 
various span stations, Moo = 2, a = 2° 
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CHAPTER 6 


CONCLUSIONS 


6.1 Summary 

Through advances in computer hardware and numerical procedures, researchers in 
the field of CFD are now able to attempt the calculation of flow about geometrically 
complex bodies. With this capability comes the need for alternatives to the single-grid 
approach to mesh generation. The creation of a single grid about a complete configuration 
or a multiple-body geometry is a time-consuming and perhaps even impossible task given 
the necessary constraints on grid quality. Researchers are exploring the use of composite 
grids for these complex applications whereby subregional grids are patched or overlaid to 
form a global mesh. The grid in each region or zone can then be generated somewhat 
independently of grids in other regions using existing grid-generation methods. 

The focus of this work is on the development of a grid-patching technique for three- 
dimensional problems that maintains global conservation. The technique is derived within 
the framework of an explicit, first-order, upwind method and is applicable in generalized 
coordinates for arbitrary point distributions on a planar zonal surface. The zonal interface 
scheme guarantees the conservation of mass, momentum, and energy across the patch 
boundary in a discrete sense. 

Calculations of the flow about a double-airfoil configuration (two- dimensional) and 
a wing-canard combination (three-dimensional) demonstrate the usefulness of the method 
in simplifying grid generation for multiple-body problems. The computations presented 
here are the first conservative zonal-grid calculations in three dimensions for patching grids 
with arbitrary point distributions at the interface. Stable, conservative calculations are 
possible, despite discontinuities in grid lines across the zonal plane, given some attention 
to the construction of the subregional grids. In the design of a patched grid system, it is 
important to ensure modest grid-density ratios across the boundary (especially in regions 
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of high solution gradient) and to choose the zone of integration such that a coarse grid 
flux results from the summation of fine grid fluxes. 

6.2 Recommendations for Further Study 

The complexity of the zonal-boundary scheme proposed in this work arises from the 
constraint that the scheme preserve global conservation at the patch boundaries. Any num- 
ber of simpler (but nonconservative) interpolation techniques could be used at patched grid 
interfaces. However, for applications in which discontinuities propagate from zone to lone 
or pass through a zonal boundary, a conservative treatment of these interfaces is necessary 
to guarantee a correct (if not smooth) solution (Benek et al., 1983; Hessenius and Pul- 
liam, 1982; Rai et al., 1984). Even in regions of smooth flow, the use of a nonconservative 
boundary technique is the suspected cause of “glitches” at a zonal interface (Rai, 1985c). 
There is clearly a need for the development of methods for the conservative treatment of 
zonal interfaces for use with composite-grid systems (both patched and overlapping). 

In theory, it is quite simple to extend the present work to the conservative patching 
of grids at an arbitrary, nonplanar surface. However, as was demonstrated in this study, 
the conservation principle is theoretically simple, but complicated to implement. There 
should be modifications of the procedure described by figure 20 to account for the case 
when zone-1 triangles and zone-2 triangles are not coplanar. This modification would 
entail projecting the zone-1 triangle vertices onto the plane defined by a zone-2 triangle 
before performing the clipping procedure to find the area of overlap. The finer the zonal 
surface is “triangulated” , the greater the degree of accuracy and computational expense. 

Another unresolved problem in three-dimensional, conservative grid patching is the 
reconcilation of a zone-1 versus a zone-2 view of a curved boundary. Consider the sketch of 
a flux-balance plane in which one edge is curved (fig. 35). This situation would occur when 
a zonal boundary intersects a body boundary. A coarse, zone-1 grid is shown in figure 35b 
(denoted by x grid points) as well as a fine, zone-2 representation of the same curved edge 
(denoted by (•) points). If the flux from the fine grid is summed to define the flux into 
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a coarse grid cell on this surface, as is recommended, then discrepancies occur near this 
curved edge. The coarse grid is not completely contained within the fine grid, and therefore 
it is difficult to make the computation conservative on the coarse grid in this region. The 
best approach to this problem is not apparent, but strict conservation in this region may 
be an unnecessary requirement. Perhaps it is not so important to be conservative near 
these edges of the zonal surface since most, if not all, boundary conditions on these edges 
are nonconservative anyway. 

The majority of the computational expense in calculating the flux interpolation coef- 
ficients is in the determination of the overlapping area between two triangles. This task 
is now done by a clipping routine borrowed from the field of computer graphics. It seems 
entirely feasible that an analytical formula could be developed to find the polygon of over- 
lap given the vertices of two triangles. An unsuccessful attempt was made to derive such 
a formula in the course of this study. The discovery of such an expression is necessary to 
permit the efficient computation of unsteady, moving grid applications (three-dimensional) 
using a patched grid system. (Moving grid calculations would require the evaluation of 
the flux interpolation coefficients at each time step.) 

The problems cited above in extending the patched grid scheme for general usage are 
all geometry related. However, for practical applications, the interfacing scheme must be 
used with second-order accurate, implicit numerical methods. The necessary modifications 
in three dimensions would directly follow the work of Rai (1985) with his two-dimensional, 
patched grid scheme. Figures 36a and 36b present first-order and second-order accurate 
(Osher method) results of Rai for the zonal blunt-body problem previously described in 
Chapter 3. Pressure contours are shown for the Mach 2 flow over a portion of a cylinder, 
and the zonal interface in denoted by the curve AB. As with the first-order results of 
this effort and of figure 36a, the second-order accurate solution is smooth and continuous 
across the zonal boundary (fig. 36b). Furthermore, the second-order accurate shock loca- 
tion agrees well with the highly accurate result of Lyubimov and Rusanov (1973) (shown 
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as square symbols in this figure). Because of the success of this two-dimensional work, 
no fundamental problems are anticipated in the use of the three-dimensional interfacing 
scheme with a second-order accurate numerical method. 

As stated previously, the zonal-boundary condition developed in this study as well as 
its two-dimensional counterpart (Rai, 1986) produce a locally first-order accurate result at 
the interface. Various analyses, including the recent work of Allmaras and Baron (1986), 
show that a given integration scheme which is nth-order accurate and used with a boundary 
scheme of n — 1-order accuracy, global accuracy still remains of order n. For this reason, 
the results of Rai in figure 36b are globally second-order accurate. At this time, however, it 
is unclear if it is possible to develop a zonal-boundary scheme that is better than first-order 
accurate at the interface. Allmaras and Baron (1986) report on a two-dimensional analysis 
of patched grid schemes in a finite volume framework using Runge-Kutta methods. The 
authors conclude that it is impossible to develop a conservative interface formulation that 
results in local second-order accuracy (first-order accuracy is possible at best). Although 
it is difficult to make a direct analogy between the schemes tested in by Allmaras and 
Baron (1986) and the one presented here, further study is warranted to determine if there 
indeed is a fundamental limitation on the accuracy of a conservative, interfacing method 
for the patching of grids. 
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(b) 



Fig. 35 View of a flux-balance plane with a curved edge 

(a) actual surface geometry (b) superimposed fine and coarse 
grid discretization 
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Fig. 36 Comparison of Rai’s (a) first-order and (b) second-order accurate 
results for a two-zone blunt body computation (Ref. 32) 
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APPENDIX A 


INTERPOLATION OVER TRIANGLES 

The following interpolation procedure is commonly used in finite-element calculations 
and is used here in the three-dimensional zonal-boundary scheme: 

1) for extrapolating grid lines into a neighboring zone to the flux-balance plane, and 

2) for interpolating the dependent variables for points on the zonal boundary in the 
neighboring zone. 

This appendix outlines the interpolation method of Chung (1978). 

We wish to interpolate for the variable u at the point (x,y) given knowledge of this 
function at the vertices of a triangle that contains (x,y) (see fig. 37). Assuming a linear 
variation of u in both x and y within the triangle, 


u = s 0 + $ix + s 2 y 


(A.l) 


where 


1 V 

x = x - 3 }^ Xi 

i=i 

1 V 

y = y- 3 


t=i 


We may determine the constants s 0 > -s i, and s 2 in equation (A.l) by constructing and 
solving the following system of equations using known values of u at the triangle vertices. 


«1 = So + SiXi + S 2 yi 
u 2 — s 0 + s l%2 + S 2 y 2 
U Z — 1 * * * 5 0 + + s 2yz 

Solving Eq. (A. 2) for s 0 > «i, and s 2 and substituting into Eq. (A.l), we obtain 


(A.2) 


U = $xUi + $ 2 U 2 + $3U 3 


M-3) 
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where 


= a» + biX + Ciy 


ai = 

\D\ {x * y3 

- xzy?) 

a 2 = 

1 

= 1^1 

(x 3 yi - xi y 3 ) 

a 3 = 


- X 

*1 = 

\D\^ 2 ~ 

yz) 

b 2 = 

1 

= \D\ 

(V3 ' 

-yi) 


f>3 = 

p^ 1 _ 

y 2 ) 

Cl = 

1 

CO 

rH |^_ 

* 2 ) 

C 2 = 

1 

= i^| 

(xj - 

-x 3 ) 


C3 = 


Xi) 

and 
















/ 1 

Xi 

yi\ 







\D\~- 

= det 

» 

X 2 

V2 









Vi 

X 3 

yzj 





It may be shown that a i = 02 = a 3 = 

To interpolate values of the dependent variables on the zonal-boundary plane, the 
variable u above is replaced by the vector Q (eq. (4.1)). 

To locate points on the flux-balance plane using this interpolation procedure, we 
must first compute the $ coefficients using the known point relationships on the zonal- 
boundary plane. Coordinate locations on the flux-balance plane are then found by the 
linear combination of these $ coefficients with coordinate values of the same triangle 
formed on the flux-balance plane. 
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APPENDIX B 


A MODIFIED SUTHERLAND-HODGMAN CLIPPING ALGORITHM 

A modification of the Sutherland-Hodgman clipping algorithm (Newman and Sproull, 
1979) is used to find the vertices of the polygon formed by the overlap of two coplanar 
triangles. A flow chart of the modified procedure is presented in figure 38, and the necessary 
changes are discussed herein. A detailed description of the original algorithm is given by 
Newman and Sproull (1979). 

We wish to find the area of overlap of the zone-1 triangle (vertices A \ , A 2, and A3) 
with the zone-2 triangle (vertices Bi, B2, and B3). (For convenience, a fourth vertex is 
added that is equivalent to the first vertex.) The area of overlap may be easily determined 
given the vertices of the polygon of overlap. The procedure for finding these vertices (Fig. 
38) involves clipping the zone-2 triangle against the edges of the zone-1 triangle. In general, 
however, it is necessary to repeat the procedure to identify all vertices, reversing the roles 
of the zone-1 and zone-2 triangles for each pair of triangles. The results of each pass are 
combined and the vertices ordered to properly describe the polygon of overlap. The area 
may then be computed using equation (4.10). 

The procedure outlined is straightforward and easily coded. The unknown vertices 
are either the intersection points of triangle line segments or are triangle vertices internal 
to the “frame” triangle (points on the visible side of a frame triangle edge). 
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PROCEDURE FOR CLIPPING AGAINST ZONE 1 TRIANGLE 



Fig. 38 Schematic of procedure: modified Sutherland-Hodgman 
clipping algorithm 
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APPENDIX C 


DETERMINATION OF RIEMANN INVARIANTS 
FOR THE 3-D OSHER ALGORITHM 


As shown in Chapter 4, it is necessary to redefine the Riemann invariants, v and u>, 
for Osher-algorithm computations on an arbitrary grid. The constraint equations (4.16) 
suggest that we choose the coefficients, c,, such that the vectors [ci,c 2 ,c 3 ] and [c 4 ,C5,c 6 ] 
are orthogonal to the vector \k x ,k y ,k z \ and to each other where 


V = CiU + c 2 v + c 3 w 

W = C4U + C5V + CqW 


(C.l) 


This may be accomplished in the following manner. 

Find a transformation matrix, T, that rotates a unit vector on the x-coordinate axis 
to the vector [k x / D,k y / D,k z / D] where D — (k x + k y + k z ) l/ 2 

jl\ l KID 
T 0 = k y /D 

\o) \K/D 

where 

( C a Cp -CpS a ~Sp\ 

T = I S a C a 0 (C.2) 

V C a Sp —S a Sp Cp J 

S a -- srn(a) = k y /D 

C a = cos(a) = {k\ + k 2 z ) 1 ^ 2 /D 

Sp = sin((3) = k z /(kl + k 2 z ) 1/2 

Cp = cos(0 ) = k x /(k 2 x + k 2 z ) 1/ 2 

and a and /? are angles described in figure 39. Then, by definition (refer again to fig. 22) 


u-kt/D = T 
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and we choose 


v = T 
w = T 




Therefore, 



and the coefficients c,- are the following components of the T matrix : 


Cj C pS^ 

C2 = C a 


C3 — S& Sp 

(C. 3 ) 

C4 = —Sp 
c 5 = 0 
C6 = C/3 

Note that if both k x and k z are zero, certain entries of the T matrix are then undefined. 
This condition is easily remedied by initially choosing to rotate the unit axis vector (x, y, 
or z) “closest” to the given {k x /D,ky/D,k z /D] vector. This practice will ensure a non-zero 
denominator in the expressions for Sp and Cp. 


93 



Fig. 39 Coordinate system and angle definitions for the determination 
of Riemann invariants 

i 
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